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Abstract Detection of 7-ray emissions from a class of active galactic nuclei (viz blazars), 
has been one of the important findings from the Compton Gamma Ray Observatory (CGRO). 
However, their 7-ray luminosity function has not been well determined. Few attempts have 
been made in earlier works, where BL Lacs and Flat Spectrum Radio Quasars (FSRQs) have 
been considered as a single source class. In this paper we investigated the evolution and 
7-ray luminosity function of FSRQs and BL Lacs separately. Our investigation indicates no 
evolution for BL Lacs, however FSRQs show significant evolution. Pure luminosity evolution 
is assumed for FSRQs and exponential and power law evolution models are examined. Due 
to the small number of sources, the low luminosity end index of the luminosity function for 
FSRQs is constrained with an upper limit. BL Lac luminosity function shows no signature of 
break. As a consistency check, the model source distributions derived from these luminosity 
functions show no significant departure from the observed source distributions. 
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1 INTRODUCTION 

Gamma-ray astronomical studies received a substantial boost after the launch of the Compton Gamma-Ray 
Observatory (CGRO) in 1991 (Kanbach et al. [19881 . Though the recently launched 7-ray missions AGILE 
& FERMI Gamma-ray Space Telescope (FGST) are expected to dramatically increase the number of 7-ray 
sources and identify new source classes, at present the 3rd EGRET (3EG) point source catalog (Hartman 
et al. 1999) provides the most complete list of GeV 7-ray sources. It contains 271 sources (>100 MeV), 
which include five pulsars, one probable radio galaxy (Cen A), 66 high confidence identifications of a sub- 
class of active galactic nuclei called blazars and one external normal galaxy, the Large Magellanic Cloud 
(LMC). In addition, 27 lower confidence potential blazar identifications are listed. Applying a different ap- 
proach, Mattox, Hartman & Reimer (1200 11 1 found 46 high confidence blazars (45 of these are present in the 
3EG catalog) and 37 plausible candidates. Sowards-Emmerd, Romani & Michelson (2003) and Sowards- 
Emmerd et al. (2004) introduced a new technique to identify 7-ray sources. Unlike earlier work (Hartman 
et al. 1999; Mattox, Hartman & Reimer 2001), where selection has largely proceeded by correlation with an 
existing radio survey, Sowards-Emmerd, Romani & Michelson (2003) and Sowards-Emmerd et al. (2004) 
attempted to obtain a more complete census of plausible blazar counterparts, sifting sources with extant 
radio survey data and then conducting a multiwavelength follow-up. Their technique ensured that even the 
plausible candidates have more than 80% good identifications. They reported 113 blazar IDs for the 3EG 
sources. 

While studying any source population, it is very important to investigate their density and luminosity 
distribution and also their time evolution. Luminosity function of a source class gives a quantitative picture 
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of the luminosity and density distribution of that source class. Blazars form the largest source class of 
identified EGRET sources, hence several papers have discussed the expected contribution of blazar emission 
to the diffuse 7-ray background. 

In order to calculate their contribution to the 7-ray background, one needs to derive the source lumi- 
nosity function. Luminosity function {(f>(L, z)) is defined as the number of sources, per unit luminosity bin 
(dL), per unit of comoving volume (dV), 

Here, we present a study of the 7-ray luminosity function of EGRET detected blazars and their evolution. 



2 LUMINOSITY FUNCTION CONSTRUCTION 



Luminosity function is constructed by binning sources in the luminosity and redshift plane from a complete 
source catalog, devoid of any selection bias. However, when the source catalog contains only a limited 
number of sources, the luminosity function can be constructed as follows. If the luminosity of a source in 
one waveband is linearly related to the luminosity in some other wavebands, one can replace the luminosity 
function of a given source class in one waveband by a scaled luminosity function from another waveband. 
This approach has been adopted in some earlier works, including Stecker, Salamon & Malkan (1993) and 
Stecker & Salamon d 19961 ), who considered linear correlations between radio luminosity and 7-ray lu- 
minosity of blazars. They approximated the 7-ray luminosity functions by the radio luminosity function. 
A similar approach has been taken by Narumoto & Totani (2006 2007) who derived the 7-ray luminos- 
ity function from scaling the X-ray luminosity function. Alternately, Chiang et al. ([1995) and Chiang & 
Mukherjee ( 1998 ) tried to construct the luminosity function directly from the 7-ray source catalog, without 
considering any correlation between radio & 7-ray luminosity of blazars. Here we adopt a similar approach 
to find the luminosity functions of Flat Spectrum Radio Quasars (FSRQs) and BL Lacs separately from the 
latest EGRET detected blazars list. We consider a hCDM cosmology ($!m=0.3, ilA=0.7). The value of 
Hubble constant (Ho) is considered to be 70 km s _1 Mpc _1 . Fig.Q]shows the overall scheme we followed 
in order to derive the luminosity function, as elaborated in the following sections. 



3 GAMMA-RAY BLAZAR CATALOG 

Chiang & Mukherjee (1998) considered the subset of 3EG catalog sources that have a minimum radio flux 
of Uy @ 5 GHz. There were 34 blazars in their list. After the publication of the 3EG catalog, Sowards- 
Emmerd, Romani & Michelson (2003 ) and Sowards-Emmerd et al. (2004 ) reported new identifications of 
EGRET sources. They calculated the over density of sources near high galactic latitudes > 20°), which 
showed a considerable number of sources when the minimum radio flux is decreased from 1 Jy to 100 mJy. 
We used the list of blazar sources from Sowards-Emmerd, Romani & Michelson (2003 ) and Sowards- 
Emmerd et al. (2004) catalog and separated them into BL Lacs and FSRQs. While constructing a gamma- 
ray blazar sample, the cutoff in the source detection significance limit has been taken as Acr for sources 
above the Galactic plane > 10°) and as 5a for sources in the Galactic plane (—10° < b < 10°). The 
final list includes 46 FSRQs and 15 BL Lacs in our sample. Three BL Lacs do not have redshift information. 
These three sources are only used to calculate the average spectral index and normalization of the luminosity 
function. 

We adopted the v v test (Avni & Bahcall [19801 ) in order to test for any source evolution. Here V 
is the volume enclosed at the known distance of the source. There exists a maximum (z max ) distance at 
which the source is at the limiting flux of the survey. V max is the volume corresponding to this maximum 
distance. Beyond z max the source flux falls below the flux limit of the survey and hence cannot be detected. 
For a system with no evolution, v v values should be uniformly distributed between and 1 and hence, 

< v v > = 0.5. Since the sample is flux limited both in radio and 7-rays, each source is assigned a z rnax = 
mm(z maXtra dio, z maxn ). While calculating z max , ra dio for FSRQs, the radio evolution function of Dunlop 
and Peacock ( fT990b is used ( f D p(z) = 10( QZ+bz2 ), where a = 1.18 & b = -0.28 ). For BL Lacs, the radio 
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Fig. 1 Flowchart of Luminosity Function Construction 



evolution function from Stickel et al. ( 119911 ) is used ( exp[^r^-}, where T(z) is the look-back time and tr 
is the evolutionary time scale in the units of Hubble time ). 

Our sample contains 12 BL Lac sources. For BL Lacs, the value of < v v > is 0.59 ±0.08. Error 

in < v v > is estimated using a = (12N)~i, where N is the number of sources in the sample. Hence 

we conclude from the y v test that BL Lacs show no measurable evolution. We also compare the y v 
distribution with a uniform distribution by KS-test and Quantile-Quantile plot, which also shows that the 
distribution is uniform. There are 46 FSRQs in our sample. For FSRQs, the value of < v v > is 0.71 ± 

0.04, which indicates strong evolution (5.2a). The v v distribution of BL Lacs and FSRQs are shown in 
Fig. |2] and Fig.[3]respectively. 

Any source class can exhibit pure luminosity evolution where only the source luminosity changes with 
redshift, while the source density remains constant. Alternately, it can exhibit pure density evolution where 
only the source density varies with redshift. More realistically, both luminosity and density evolution can be 
expected. Considering the limited number of 7-ray sources, it is difficult to examine an evolution model that 
incorporates both luminosity & density evolution. We examined two pure density evolution models, [(1 + 
zf 1 and exp(/3 2 T(z))]. We found large errors in the density parameter values as derived from our source 
list (Pi = 5.8^15 and — 10.9^2'g)- These large errors prevents us from drawing useful conclusions. 
Since pure luminosity evolution is more often observed at other wavelengths, we examine such a model for 
7-ray blazars. 
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4 LUMINOSITY EVOLUTION FUNCTION FOR FSRQS 

v analysis shows a clear indication of evolution of FSRQs in 7-rays. We consider the pure luminosity 
evolution (density of sources is constant with z) of these sources. The luminosity of a source at a redshift z 
can be written as 

L{z) = L x f(z) 

where Lq is the luminosity at zero redshift and f(z) is the luminosity evolution function. 



(2) 



Two types of luminosity evolution functions, exp(— ^) and (1 + z)' are considered. Here T{z) denotes 

- method to find the evolution parameters r and j3. For the 



look-back time. We used the modified 

optimum parameter value, v v should be uniformly distributed between and 1, thus < v v >= 0.5. 
For very high values of the evolution parameter, we sometimes get z m i n as well as z max . The physical 
significance of z m i n is that, for that particular evolution parameter value, the source cannot be observed 
below z m i n . For all cases with z max > 5 we assumed z max — 5. Fig.Blshows the variation of < , r v > 



with different r values. Horizontal dashed lines show the 1 a error in < y^- 



>. For each value of the 



evolution parameter (r), the distribution of v is compared with a uniform distribution. We find t = 

0.16 ± 0.02. KS-test is performed which shows the distribution of v v is uniform for r = 0.16 ± 0.02. 

We randomly took 46 points from a uniform distribution and studied the Quantile-Quantile plot of the v v 
(after de-evolution for r = 0.16) with these 46 randomly chosen points. It shows no significant departure 
from linearity and thus again demonstrates that the v v distribution is uniform. Fig.[3]shows the variation 

of < jr— > with f3 considering the second evolutionary model (power law). We find [3 



which < 



3.0±° j for 



> becomes 0.50 (within 1 a). KS-test & Quantile-Quantile plot also show the distribution 



is uniform for (3 = 3.0^0 1- 

For both models of evolution, all luminosities have been de-evolved to z = and these are used to 
construct the de-evolved luminosity function by the standard nonparametric — method and the likelihood 
method. We used the average spectral index of the sources for luminosity function determination. 

5 SPECTRAL INDEX DISTRIBUTION OF FSRQS AND BL LACS 

The photon spectral index distribution with redshift of FSRQs and BL Lacs from our source list are shown 
in Fig. [6] The photon spectral index distribution of these objects with luminosity are shown in Fig. [7] In 
Fig- El no evolution of FSRQs and BL Lacs have been considered. Of the 15 BL Lac sources, the three BL 
Lac sources without redshift information are considered to have the average redshift of the distribution. To 
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Fig. 6 Spectral index distribution with redshift. The 
circles are the FSRQ objects and the diamonds are 
the BL Lac objects. 



3.5 




1 ' ' 1 

44 46 48 < 50 

log 10 (Luminosity (erg s )) 

Fig. 7 Spectral index distribution with luminosity. 
The circles are the FSRQ objects and the diamonds 
are the BL Lac objects. 



derive the average 7-ray spectral index, we assume that the intrinsic spectral index distribution (ISID) can 
be described by a Gaussian (Venters & Pavlidou (2007]), 

lSlD(a)da = — exp 

V27TCTO 

The likelihood function of spectral index distribution has been adopted from Venters & Pavlidou (2007 ). 

Fig.[8]and Fig.|9]show the 68%, 95% and 99% likelihood function contours of ao and do f° r FSRQs and 
BL Lacs respectively. Venters & Pavlidou (2007) also calculated the average spectral index for FSRQs and 
BL Lacs. However, their source selection criteria did not require > 4a detection of sources when averaged 
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Fig. 8 68%, 95% and 99% likelihood function con- 
tours for FSRQs 



Fig. 9 68%, 95% and 99% likelihood function con- 
tours for BL Lacs 



Table 1 The Average 7-ray Spectral Index Value for Blazars 





This work 


Venters & Pavlidou (2007) 


FSRQ 


2.34 ±0.15 


2.30 ±0.19 


BL Lac 


2.08 ±0.18 


2.15 ±0.28 



over the first four phases of EGRET observations. The values of «o & &o of our calculations and from 
Venters & Pavlidou ( 20071 are given in Table [JJ 

6 LUMINOSITY FUNCTION CONSTRUCTION OF FSRQS 

0, we constructed their luminosity function by the 



After de-evolution of the source luminosities at z 
following methods. 

(a) The luminosity function can be written as 



1 N 

i=l 



dL ^ V max (i) 



(4) 



Here, Lq is the de-evolved luminosity. For our data set, 



method only gives the upper-end index of the 



luminosity function. The break luminosity (if any) and the lower end index cannot be found by this method, 
(b) Maximum Likelihood function has been generated for the distribution of sources, considering a broken 
power law luminosity function. The high luminosity end slope has been taken from the non-parametric 
— method (as discussed in (a) ). The break luminosity and the lower end index are the free parameters, 
and they are constrained using the likelihood function. 



6.1 



Method 



The non-parametric — method has been employed in order to find FSRQ luminosity function (cf>(Lo)). 
FSRQs have been binned into different luminosity intervals. The high end part of the luminosity function is 
well fitted by a power law with index of 2.5± 0.2 (assuming exponential evolution function with r = 0.16). 
For power law evolution function ( j3 — 3.0), the best fitted value of the power law index is 2.4 ± 0.2. At the 
low luminosity end there is a hint of a turn over, but data points are insufficient to confirm it unambiguously. 
Fig. [TO] & Fig. [TJJ show the variation of 4>(Lo) (determined by — method) with de-evolved luminosity 
considering the exponential evolution function and power law evolution function respectively. 
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Fig. 10 De-evolved luminosity function (</>(Lq)) of 



FSRQs (r = 0.16) (, 
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Fig. 11 De-evolved luminosity function (</>(Lrj)) of 
FSRQs (fi = 3.0) method). <j>(L ) is given in 

(^r) x {^-) 3 per luminosity (10 45 erg s _1 unit) in- 
terval. 



In order to examine a possible change of index at the low luminosity end and to find the break luminosity 
(if any), we return to the source redshift distribution. Maximum likelihood analysis is used to derive the 
break luminosity and the low-luminosity end index. 

6.2 Maximum Likelihood Analysis 

We assumed a broken power law form of the de-evolved luminosity function. 



Ln < Lf 




(5) 



Here 4>(Lq) (= d y^ Lo ) is the de-evolved luminosity function and cf>o is the normalization of the luminosity 
function. Lb is the break luminosity. We fixed the high luminosity end index (o^) of the luminosity function 
from — method (as described in the last subsection). The redshift distribution of EGRET detected FSRQs 
is used to find both break luminosity {Lb) and the power law index {ax) of the low-luminosity end of the 
4>{Lq). The likelihood function for the redshift distribution of FSRQs (similar to Chiang & Mukherjee 1998 ) 
is given by 



N 



[dN/dz]i(zi 



J* ma *[dN/dz]i(z)dz 



Source distribution with z is given by, 

~dN 



dz 



(z) = @(z) 



dV 

dz 



4>{L )dL 



Here L i<Hm (z) is given by, 



{i + z y-^f{z) 



(6) 



(7) 



(8) 
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f{z) is the 7-ray evolution function and a 7 is the average 7-ray energy spectral index of FSRQs. Fi t u m 
is the limiting flux of the ith source (Chiang et al. 119951 1. &(z) is the radio completeness function as defined 
in Chiang & Mukherjee ( 1998 ), which is the fraction of the FSRQs at a given redshift z having radio flux 
greater than the radio limiting flux of the survey (100 mJy). 



Pi,,. 



dN ( f°° dN 
dP—{ / dP— 



(9) 



where, — cx [(^-)° 83 + (is-) 1 ' 96 ] 1 [Dunlop and Peacock (1990)] 
dP P Pb Pb 

HereP B = 10 25 26 WHz^sr 1 



i.lim.radio 



Flim,radioD L (z) 



(i 



io fDp(z) 



(10) 



The de-evolved limiting luminosity Pi,u m ,radio(z) is a function of z and varies from source to source. The 
radio evolution function is taken from Dunlop and Peacock dl9901 l. 

For an exponential evolution function, we find a break luminosity (Lb) of 1.7 x 10 46 erg s . This 
analysis gives an upper limit to the low luminosity end power law index ot\. Fig. [T2] gives the 99% upper 
limit of a\ as 1.3 and the 95% upper limit as 1.1. For the power law evolution function, we find break 
luminosity Lb is 5.8 x 10 46 ergs _1 . Fig.[l3]gives the 99% upper limit of ct\ as 1.1 and the 95% upper limit 
as 0.9. 
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Fig. 12 68%,95%,99% contours considering ex- 
ponential evolution function. Dashed line shows the 
maximum likelihood L b value 



Fig. 13 68%, 95%, 99% contours considering power 
law evolution function. Dashed line shows the maxi- 
mum likelihood Lb value 



The derived luminosity function parameter values are given in Table|2] The source distribution (dN/ dz) 
has been over plotted with the model distribution for both exponential and power law evolution model in 
Fig. [14] For the model luminosity functions the 95% upper limit of low luminosity end index is taken. 



7 LUMINOSITY FUNCTION FOR BL LACS 



< 



> study shows no evidence of evolution for BL Lacs (Fig. [2]). We studied the variation of < 



> 



V mfx ' " * «'>"">™>"«»-«'«"'«»»'»»"»"-"»"V'W — ^ V, 

with the evolution parameter (t) of a pure luminosity evolution function (similar to FSRQs). It is found that 
< v v > is not sensitive to the r variation (Fig. [15). We construct the luminosity function of BL Lacs by 
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Fig. 14 The histogram is constructed from our FSRQ source list. The thick solid line shows the 
distribution for the exponential model. The thick dashed line shows the distribution for the power 
law model. 



Table 2 Luminosity Function Parameters 





Exponential evolution function (t) 


Power Law evolution function (J3) 


L B (in 10 46 ergs" 1 ) 

ai (95% upper limit) 
qi (99% upper limit) 


1 7+ 1 ' 8 

l ''-0.A 

2.5±0.2 
1.1 
1.3 


r o+2.5 
->-»-0.9 

2.4±0.2 
0.9 
1.1 



following similar methods used in the case of FSRQs. We find a single power law luminosity function with 
an index of (2.37 ± 0.03) by method (Fig. [Hi. 

Since the source list of BL Lacs is too small to independently determine the lower end luminosity index 
and break luminosity for BL Lacs, we assumed values similar to those derived for FSRQs. 

7.1 Average EGRET Limiting Flux 

The EGRET flux limit was not found to be identical for all sources. The empirical relationship between flux 
limit and statistical significance is given by (Chiang et al. |1995l ). 
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Fig. 15 Distribution of < y±— > values for BL Fig. 16 <j>(L) vs L plot for BL Lacs (,7^— method). 
Lacs with different r values. Dashed horizontal lines i s g i ve n in ) x (-2a) 3 per luminosity (10 45 
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indicate the 1 a error in < 



v 



While calculating the V max of each source and also while finding the low end luminosity index and break 
luminosity, the appropriate flux limit for each source has been calculated. In order to compare the observed 
density distribution with the model distribution, the average limiting flux (~ lx 10 -7 phcm -1 s _1 ) is used. 

8 NORMALIZATION CONSTANT OF LUMINOSITY FUNCTION 

The normalization of the luminosity function has been found by integrating the luminosity function over 
the redshift range zero and z max , and above the EGRET limiting luminosity. 

N obs = / — dz / <f>(L )dL (12) 

J0 dz JL lim (z) 

Here, N f, s is the number of FSRQs (BL Lacs) detected by EGRET. Since each EGRET source has a differ- 
ent limiting flux, in order to calculate the limiting luminosity Lu m (z), the average limiting Qax(Fn mtave ) 
has been used (1 x 10~ 7 phcm~ 2 s _1 ). The limiting luminosity is given by 

4tt x D 2 L F Um , ave 
(1 + z) 2 a x f(z) 

Here a is the average photon spectral index of EGRET detected FSRQs (BL Lacs) and f(z) is the 7-ray 
evolution function of FSRQs. For BL Lacs, f(z) is unity. Also since the coverage in the southern hemisphere 
is 0° < b < —40°, a correction factor for the loss of solid angle is included. 

9 DISCUSSION 

We investigated for the first time the luminosity function and evolutionary nature of FSRQs and BL Lacs 
separately in 7-rays. The construction of the luminosity function is of great importance in order to under- 
stand the nature of these populations as a whole. It also plays an important role in estimating their con- 
tribution to the 7-ray background. We find no evolution for BL Lacs, whereas for FSRQs, there is strong 
indication of evolution in 7-rays. Padovani et al. 2007 found the luminosity function of BL Lacs and FSRQs 
using radio and X-ray data. They also found no evidence of evolution for BL Lacs, whereas FSRQs show 
strong evolution. For FSRQs, we assumed pure luminosity evolution. Two types of evolution functions 
(power law and exponential) are considered in this work. The source luminosities are de-evolved to zero 
redshift and the de-evolved luminosity function is constructed by y± — method. This analysis only gives 
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the high luminosity end power law index of the luminosity function. Assuming a broken power law model, 
we used likelihood analysis of source density distribution to find the break luminosity and low luminosity 
end index. Only an upper limit for the low end index is derived. For BL Lacs, the luminosity function is 
described by a single power law with inadequate data to derive a break luminosity. 

In earlier work, Chiang & Mukherjee ( 1998 ) found a much flatter high end luminosity index (a.2=2.2) 
than our estimation (a 2 =2.5 for exponential model and 2.4 for power law model). They also found a lower 
break luminosity (Lb = 1.1 x 10 46 erg s^ 1 ) than our estimation (given in Table[2]). Contrary to this work, 
they studied the FSRQs and BL Lacs as a single source class. We constructed the luminosity function and 
evolution of FSRQs and BL Lacs with a source list having almost twice the number of sources than Chiang 
& Mukherjee ( 1998). We used a AC DM cosmology, while they considered an Einstein-deSitter universe. 
Though we do not have any quantitative measure to select between the two forms of the evolution function, 
the source distribution (Fig. n~4b suggests that the exponential evolution function is a better representation 
than the power law evolution function. 

The low end luminosity function is not well constrained due to a lack of adequate number of sources. 
Considering the much higher sensitivity of FGST, one expects a clearer picture to emerge of the distri- 
bution of these source classes from FGST data. The blazar contribution to the Extragalactic Gamma-Ray 
Background (EGRB: see, Sreekumar et al. 1998, Strong, Moskalenko & Reimer 2004) together with the 
contribution from AGNs with different inclination angles, will be presented elsewhere. 
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